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, _^ ■ Abstract 

^N ■ This paper studies the numerical solution of traveling singular sources problems. In such 

problems, a big challenge is the sources move with different speeds, which are described by 

.^ ' some ordinary differential equations. A predictor-corrector algorithm is presented to simulate 

the position of singular sources. Then a moving mesh method in conjunction with domain 
decomposition is derived for the underlying PDE. According to the positions of the sources, 

i-Q , the whole domain is splitted into several subdomains, where moving mesh equations are solved 

C^ ' respectively. On the resulting mesh, the computation of jump [u] is avoided and the discretization 

of the underlying PDE is reduced into only two cases. In addition, the new method has a desired 
second-order of the spatial convergence. Numerical examples are presented to illustrate the 
convergence rates and the efficiency of the method. Blow-up phenomenon is also investigated for 

S^ ' various motions of the sources. 

fvj ■ Keywords: Moving mesh method; Domain decomposition; Traveling singular sources 

00 " 

^ ■ 1 Introduction 
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CN ' We take the one-dimensional moving singular sources equation 
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ut — Uxx = 2_^Fi{t,x,u)5{x — ai{t)), — oo < a; < oo, t > 0, (1) 

i=0 

C^_' u{x,0) = uo{x), — oo < X < oo, (2) 

u{x,t) — )■ as |x| — )■ oo, t > 0. (3) 

as the model problem in this paper. Here (7 > is the number of singular sources. The initial value 
uo{x) is taken to be continuous and compatible with the boundary conditions, i.e. uo{x) — )• as 
|x| — )• 00. The local source functions Fi{t, x,u) (i = 0, 1, . . . , g — 1) might be given a priori or can be 
determined from some additional constraints on the solution. The traveling sources are located at 
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ai{t), i = 0,1, ... ,q — l. In general, their velocities can be described by several ordinary differential 
equations 

-^ = ^lji{t,ai{t),u), i = 0,l,...,q-l, (4) 

which are coupled with the solution u. We assume that the sources do not intersect with each other 
during the time in consideration. This model arises in many areas such as laser beams traveling 
problems where u is the temperature of the material |lj, or free-boundary solidification problems 
where ai{t), i = 0,1, ... ,q — 1, are the moving interfaces between different phases [2]. 

It is well known that the solution of the model is continuous and piecewise smooth [3j. However, 
the derivative of the solution has a jump at each source due to the delta function singularity on it, 
and the jump is given by ^ 

[ux]{a,{t),t) = -Fi{t,ai{t),u{ai{t),t)), i = 0,1, . . . ,q - 1. (5) 

This leads the standard numerical methods, either finite difference method or finite element method, 
might fail when crossing the time-dependent source positions. Various approaches, have been used 
to deal with the delta function singularity, such as the immerse boundary (IB) method and the 
immerse interface method (IIM) [2l[3l|5H8]. For the IB method originally proposed by [9j, the delta 
function is approximated by an appropriately chosen discrete delta function. Beyer and LeVeque [2] 
studied various cases of the model ([I]) — ([3]) for the IB method with q = 1, and the source position 
ai{t) being priori specified. In contrast, the IIM first introduced by LeVeque and Li [B] incorporates 
the known jumps of solution or its derivatives into the finite difference scheme to obtain a modified 
discretization scheme. Li [3] developed an IIM numerical algorithm on the uniform mesh for the 
model dU)-© and Q with q = 1. 

The model (P) — ([3]) and (JH is more difficult to be solved when the source function Fi{t,x,u) 
is high nonlinearity. In this case, the solution is always blow-up in some finite time T > if 
the sources are stationary or move at sufficiently low speed, while blow-up will be avoided if the 
sources move at sufficiently high speed (see e.g. [Il[l0l[Tl] ) . As the solution evolves singularity, the 
uniform mesh method always become computationally prohibitive. Hence, moving mesh method 
has to be employed, which is one of the most popular adaptive methods and have been successfully 
used to investigate the blow-up phenomenon |121ll3j. In MMPDE's approaches, the movement of 
the mesh is controlled by the moving mesh partial differential equations (MMPDEs) based on the 
equidistribution principle [S]. Among these MMPDEs, MMPDE4, MMPDE5 and MMPDE6 are 
popular to use. Readers interested in the moving mesh method and its applications can refer to the 
books [lailS]. 

Recently, several papers have been devoted to moving mesh method for the model ([I|) — (l3|) with 
a priori specified source position ai{t) for q = 1 [Hll7j and for g > 1, in which the sources move with 
the same speed [T8l[19j. This paper is a further study of [19] and [2^ for the model ([I]) — ([3]) with 
general movement of the sources, which do not intersect with each other during time evolution. First, 
we choose a finite observed domain containing all sources with appropriate boundary conditions, 
and divide it into q + 1 subdomains by the q sources. Obviously, the sizes of the subdomains 
are changed as the sources traveling. MMPDEs are applied on each subdomain to obtain a local 
equidistributed mesh on it. Then the underlying PDE (JT]) is solved on the whole observed domain 



with the mesh composed of the local mesh on each subdomain. Taking the advantages of domain 
decomposition [21], MMPDEs could be solved efficiently by parallel computing. Moreover, It can 
be found that the computation of [u] is avoided, thus the discretization scheme for the underlying 
PDE becomes very simple. In addition, our method has an expected second-order convergence in 
space. 

The organization of the paper is as follows. In section [21 we introduce the moving mesh method 
in conjunction with domain decomposition for the model problem. In section [21 the discretization 
schemes for the physical problem will be derived in detail. In section [H several numerical examples 
are given to demonstrate the numerical efficiency and accuracy of our method. The conclusions are 
presented in the last section. 

2 Moving mesh method in conjunction with domain decomposition 

In the last decade, moving mesh method in conjunction with a Schwarz domain decomposition has 
been developed by Haynes and his co-workers (see i.e. |22|, I23) and references therein). And in this 
section, we will introduce a slight different method, that is, moving mesh method in conjunction 
with a non-overlapping domain decomposition. 

Denote the observed domain by [x/, 2;^.] and assume it containing all sources, that is, xi < ao(0 < 
ai{t) < ■ ■ ■ < ag-i{t) < Xr- Here xi, Xr are either constants or variables of t. Then, the observed 
domain is divided into q + 1 subdomains [aj_i, aj] (i = 0, 1, . . . , g) with a_i = xi, Ug = Xr, by the 
q sources respectively. Obviously, the sizes of the subdomains are variables of t too. 

Let X and ^ denote physical and computational coordinates, respectively. Without loss of gener- 
ality we assume the computational domain is [0, 1]. Then an one-to-one coordinate transformation 
between the observed domain [x;,a;j.] and the computational domain [0, 1] is defined by 

x = x{^,t), ^£[0,1], (6) 

with 

x{0,t) = Xi, x{l,t) = Xr- 

For a given uniform mesh, ^j = -^, j = 0,1, . . . ,N, on the computational domain, the corresponding 
mesh on the observed domain [x/, j;^] is 

Xi = Xo{t) < Xi{t) < • • • < X7V-l(t) < XN{t) = Xr- 

In our method, the coordinate transformation Q is determined as a piecewise smooth function. On 
each subdomain [ai_i, Oj], i = 0,1, - - - ,q, it is the solution of an MMPDE which is derived from the 
equidistribution principle. In the literature, the following MMPDEs 

1 d ( ,,dx\ 
d'^x 1 d / dx\ 



which known as MMPDE4, MMPDE5 and MMPDE6, respectively, are popularly used after they 
were originally established and analyzed in |14j . Here M = M{x,t) is the monitor function giving 
some measure of the solution error on the physical domain and r > is a parameter representing 
a timescale for adjusting the mesh toward equidistribution. In the asymptotic case t — )• oo, the 
solution of MMPDE4, MMPDE5 and MMPDE6 would satisfy the equidistribution principle, which 
is stated that [14J 






0. 



(10) 



For more details about MMPDE, one can refer to [Tl] or the recent book [TB]. In this paper, 
MMPDE6 dHD with the boundary condition 



^(Of_i'0 = cti-iit)^ x{^js,t) = ai{t), 



0,1,..., q, 



(11) 



is employed as an example to describe our moving mesh strategy in conjunction with domain de- 
composition. Here j* is some fixed index satisfying < j| < N. The resulting mesh, used to solve 
the model problem on [x/, x^], satisfies the property that a fixed mesh point is located on each source 
during the time in consideration, i.e. Xj^ = ai{t). 



Given the old mesh .t^^ on the observed domain 
[a;; , aJr] and tlie corresponding solution on the mesh. 



Compute the monitor function A/ on the mesh. 



Solve MMPDE6 © with the boundary condition ITTJ 
on each subdomain [ttz-i, cii]. i — 0,1, ... .q. 



Combining the local equidistributed mesh on each 
subdomain to give the new mesh x''^^ on [xi,Xr]. 



Figure 1: The moving mesh strategy in conjunction with domain decomposition. 

Figure [1] shows the moving mesh strategy in conjunction with domain decomposition. Here the 
computation of the monitor function will be presented in section 21 And MMPDE6 ([9]) is solved by 
the following finite difference scheme 
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in our numerical examples, where At„ = tn+i — tn and M, 



i+ 



(Mj+i + M,)/2. 



The new mesh could be obtained very efficiently by parallel computing based on domain decom- 
position methods [2IJ- And it is best in the sense of equidistribution on each subdomain. On the 
other hand, we will found in the next section that the computation of the jump [ii] is avoided, hence 
the discretization scheme for the physical PDE ([1]) becomes very simple. 



3 Model discretization and final algorithm 

In this section, we derive the discretization schemes for the physical model problem ([I]) — ([3]) and dH) 
on the observed domain [x;, x,-] with appropriate boundary conditions. Then present a full algorithm 
of moving mesh method for the model problem. 

3.1 Discretization schemes 

For an arbitrary function / = f{x,t) = f{x{S,,t),t), we have 



f = ^{x{tt),t) 



= ft + fxX. 
^ fixed 



Through the coordinate transformation ^, we can rewrite equation ([TJ on the computational 
coordinates as 

q-l 

ii - UxX -Uxx = ^^ Fi{t, X, u)6{x - ai{t)). (13) 

Since the right-hand side of (fT3|) vanishes when x / aj(t), that is, 

U — UxX — Uxx = 0, (14) 

we conduct the discretization scheme for (jlSp on the above equation as [19j, with each term on the 
left-hand side of (|14|) containing the information of jumps when they cross the sources. Physically, 
the value u of ith source changes smoothly as time evolution, which means the jump of the directional 
derivative of u{x,t) along the vector {a'-{t), 1) is zero [JlfTO]. i.e., 

[Ut](a,{t),t) + [Ux]{a,it),t)<^'M = 0, i = 0,l,...,q-l. (15) 

Recalling that Xjs = ai{t), it follows that 

[u](a,{t),t) = [Ut + UxX](^a,(t),t) = [ut](a,{t),t) + iK](a,(i),i) =0, i = 0, 1, . . . , g - 1. (16) 

By using the above equation, we can deduce from (fn|) that 

[Uxx]{a,{t),t) = [u- Ux±](a,{t),t) = N(a,(t),t) " ^Wx]{a,{t),t) = -"i(*)K](a,(i),i) , (17) 

i = 0,1, . . . ,q — 1. Then we obtain immediately 

[uxx](a,{t),t) =i}i{t,ai{t),u)Fi{t,ai{t),u{ai{t),t)), i = 0, 1, . . . , g - 1, (18) 

by taking (0]) and © into (fTTl) . 

Similarly to p^, the discretization scheme for (jlSp are divided into two cases due to Xj^ = ai{t) 
during time integration. For j / j|, i = 0, 1, . . . , g — 1, the mesh point not located at the source, p^ 



is discretized by standard center difference for spatial variable and backward difference for temporal 
variable, that is, 

where W = x'^ — x"-_^. Here x^, u^ are the mesh and the solution on it at time step tn, respectively. 
For j = jf, the mesh point just located at the source, the jump informations should be incorporated 
into the discretization scheme. For this case, the discretization scheme for (11311 reads 






2 -- ' ,n+l^ 



^^i«= ) = 0, (20) 

where V^' « V'*(in+i, <+',«"+'), i^^K^+') « F^(tn+i, ar+\ <.+'), i = 0, 1, . . . ,g - 1. 

In the above schemes, we need the source position aj(t„+i) at time step tn+i- For the general 
movement @, it is computed by the following Crank-Nicolson scheme 

«r' = «r+^(c^'+c), i = o,i,...,g-i, (21) 

as in [20j. If ipi{t,ai{t),u), i = 0,1, . . . ,q — 1, are independent of u, the source position a^^ and 
the speed ■0""'' can be calculated in advance before solving the discretization schemes for MMPDE6 
([9|) and physical PDF (|13p . Otherwise, the resulting system would be too complicated to be solved. 
In this case, we decouple the discretization system by a predictor-corrector algorithm. For the 
predictor step, assume ■0""'' = i^f and solve ([2T]) to get an approximate variable a* of a""*" . Then 
substituting t/^""'' and a* into the discretization schemes for Q and (jlSp to obtain an approximate 
solution u* of 'u"^-'^. For the corrector step, compute ip^'^ = ■ilji{tn+i,ai,u*) and solve (f2T]) to get 
a""*" , then obtain the solution u^~^^ at time step tn+i by the discretization schemes for ([9]) and (|13p. 
To complete the discretization schemes, we require an appropriate condition for u on the bound- 
ary of the observed domain [j;;, x,.]. For the observed domain is small enough, we employ a third-order 
local absorbing boundary condition (LABC) proposed in [24j 

SsqUj: + Uxt ± Sq^/s^U ± 3^/s^Ut = (22) 

for ([T|) as in |19p20j. Here sq is an user-defined parameter, the plus sign in "it" corresponds to the 
LABC at the right boundary Xr, and the minus sign corresponds to the one at the left boundary xi. 
Under the map ([6]), we get the LABC for (fT3|) as follows 



iix ± 3^/s^u + SsqUx - UxxX ± sq^/s^u ± {-Z^^Uxx) = 0, (23) 

where the plus sign in "it" corresponds to the right boundary, and the minus sign corresponds to 



the left boundary. According to |19p20j. a finite difference scheme for (|23p is 



1 / ur' -ulY u^{-u-,\ ^ ^ur-K , J , ^x^^' - xn n 



) - 3^^^^^^ + 3 (.0 + Vi^^ 



"1 



So^/5^«r' = 0' (24) 



^n+l ™n\ „n+l _ 0,,"+l i ,,"+1 



on the left boundary, and 

i / "Af+1 "AT-I "Af+l "TV-l \ , o / — ^iV ~ ^N , n I — ^N ~ ^N \ "-N+l "-N-l 



+ ^V^ A. + 3 Uo - ^/i^- 



At„V 2/1^+1 2/1^ y . ^v-u ^^^ ,^^.u v-u ^^^ y 2/,^+i 

-[ At^ ) i^;^^ +'°^"° -°- ('') 

on the right boundary. Here two ghost points x_i and xn+i are used. On the other hand, if the 
observed domain is big enough or else, Dirichlet boundary conditions are employed. 

3.2 Full algorithm 

We close this section with a full algorithm in Figure [2] for the model problem ([Ij — (l3|) and @. Here 
the choice of the time step At„ will be determined in the following concrete examples, and Tol > 
is set to be 10"^^. 

4 Numerical examples 

In this section, we present some numerical examples to verify the convergence rate and illustrate 
efficiency of the full algorithm in Figure [2j 

Example 1. We consider a nonlinear moving interface problem with the following exact solution 



sm{u}ix)e '^1*, X < ao{t), 

sin(c<j2(l - x))e~'^2*, x > ao{t), 



^(^'*) = \ . , ,. ^^ „<^2t . ,.. (26) 



for some choice of ui and UJ2- The interface ao{t) is determined by solving the scalar equation 

sin(a;iao)e"'^?* = sin(w2(l - ao))e~'^2*^ (27) 

so that u{x, t) is continuous across the interface. 

The equation (j27|) has a unique solution on [0, 1] if we take, for example, vr < wi, U2 < 2tt. Then 
we have the ordinary differential equation for the motion of the interface 

doo ^ {ujl-uj^)u{ao,t) 
dt UxiaQ,t) -Ux{a^,t) 



^ Prepare the initial values x^\ u^^ and the 




Figure 2: Full algorithm for numerical solution of the model problem ([T|-(I3|). 



Based on the jump conditions, the source function Fo(t,2;,n) is 

Fo{t,x,u) = -[ux\ao = Ux{aQ,t) -■u^(a(j',t) 

= uji cos(cjiao)e~'^i* + 0J2 cos(cj2(1 - ao))e~'^2*. 



(29) 



Same as in [2], we take wi = 57r/4, a;2 = 77r/4. The observed domain is set by [0, 1], where the initial 
position of the interface is ao(0) = 0.58333. Since we have the exact solution, Dirichlet boundary 
conditions are employed. In this example, we simply use the uniform time step, i.e. At„ = const, 
and the total number of the time meshes is L. The monitor function for MMPDE6 ([9]) takes the 



form 



M{x,t) = (1 



du 
dx 



+ 0{{x-ao{t)f + e) 



-1/4 



(30) 



where < 9 < l,0<e^l. This is consistent with the choice in |19p20j . In practice, smoothing the 
monitor function can improve the accuracy of the numerical solution, and we utilize the smoothing 
technique proposed in |25] . Here the parameters in MMPDE6 ([9]) and the monitor function (j30p are 
given by r = 10-^ 6 = 0.5, and e = lO^/A^^. 

Since backward Euler scheme is used to solve the physical PDE in this paper, the truncation error 
for time discretization is only first-order. To verify our algorithm has a second-order convergence 
rate for space, the number of L should be fourfold when A^ is double in the convergence test. 
Computational results with different number of A^ and L at the time T = 0.1 are listed in Table [TJ 
where the errors are defined as 



En,l 



U -Ue 



En,l 



U -Ue 



rpOL 



"0 - "o 



Here, Ue is the true solution, Uq used as the exact interface is the solution of a zero- finding MATLAB 
function /zero for (|27p . The numerical solution U is obtained by the algorithm where ao(i) and ao(0 
are exactly calculated. And C/, cxq represent respectively the numerical solution and interface, solved 
with the full predictor-corrector algorithm. The ratios in Table [T] are E2n,4:l/E]\[,l, E2N,iL/EN,L 
and -E'^iV4L/-^^'L' respectively. It is shown that our algorithm solves the solution and the interface 
very well, and has a second-order convergence rate for space, i.e. 0(1/A^^). Additionally, compared 
the corresponding results in [3j, our algorithm is better than the method proposed in [3]. 



Table 1: Error and convergence rates at T = 0.1. 

ijjv.i ratio Ejyj^ ratio E'^ ^ r; 



40, 40 


1.0931O-2 


- 


1.3721C-2 


- 


8.2720e-3 


- 


80, 160 


2.6996C-3 


0.24697 


3.3908e-3 


0.24713 


2.0540e-3 


0.24830 


160, 640 


6.6945C-4 


0.24798 


8.4144e-4 


0.24816 


5.1038e-4 


0.24848 


320, 2560 


1.6687e-4 


0.24927 


2.0981e-4 


0.24935 


1.2732e-4 


0.24947 


640, 10240 


4.167Se-5 


0.24976 


5.2408e-5 


0.24978 


3.1807e-5 


0.24982 


1280, 40960 


1.0416e-5 


0.24992 


1.3098e-5 


0.24993 


7.9500e-6 


0.24994 


2560, 163840 


2.6038O-6 


0.24998 


3.2743e-6 


0.24998 


1.9874e-6 


0.24998 
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Figure 3: Mesh trajectories and the profiles of u in physical variable and computational variable (from left 
to right) as time changes with N — 24. The solid lines are the computed solution and the dots are the exact 
solution on the mesh. 



Figure [3] presents the profiles of tlie solution in physical variable and computational variable and 
the evolving mesh from t = to t = 0.1. The number of the mesh is A^ = 24, with half mesh points 
on each side of the interface. We can see that we get excellent resolution of the example even with 
a grid as coarse as A^ = 24. 

The rest examples are from traveling heat sources problems with the solution may be blow- 
up [ll[T71[19l[20] . If not specifically pointed out, the initial value is given by 



u{x, 0) 



cos^(7rj;/2), 
0, 



-1 < X < 1, 

otherwise, 



(31) 



the observed domain is set by [—10,10] with n(— 10, i) = u(10,t) = 0, and the source functions 
Fi{t,x,u) are simply specified by 



Fi{t,x,u) = l + u^, 



0,1,..., g-1. 



(32) 



The resulting nonlinear system is solved by Newton iteration with the tolerance tol = 10 
The monitor function for MMPDE6 ([9]) takes the form 



Mix,t) = e,+iuP + e, 



du 



dx 



q-l 



+ Y^9,{{x-a^{t)f + e)-^/\ 



(33) 



j=0 



where the parameters < 6*4 < 1, X]?=o ^i = l!0<e^l,p>0 will be determined later. For 
non-blowup case, the following graded time steps [HfTT] 



n 



L 



n = 0,l,...,L, 



are used with [0, T] the time integration interval and L the number of time meshes. While for 
blow-up case, the time step Af„ = tn+i — tn is chosen to be 



Atn = min < 



f^, 



fJ- 



IXLBjX-i t u/j 



where e is same in the monitor function, /i is a small positive constant with fi = 10^^ in the test. 

Example 2 (Linear moving sources). We consider that all sources move with a constant velocity k, 
and the position of the i-th source has a constant distance di to the 0-th source, i.e., 



a'i{t) = k, ai{0) = di, 



0,1,. ..,5-1, 



where do = 



10 



Since our method is trivial for multi-sources case, only q = 1,2 are considered. Different velocities 
k are investigated in [11[T71[19] and we specify k = 2 here. With this velocity, blow-up would occur 
for q = 2, and be avoided for q = I. The parameters are set by r = 10~^, 6*0 = 0.9, 6i = 0.1, and 
e = lO^iV^ for g = 1, while r = 5 x 10"^, 60 = 61 = 0.3, 65 = 0.4, p = 2, e = 10-^ and di = 2.5 
for g = 2, respectively. 

The profiles of the computed solution in physical variable and computational variable and the 
evolving mesh are presented in Figure 2] for q = 1 and in Figure for q = 2. For simplicity, each 
subdomain has 50 mesh points, i.e. N = 100 for q = 1 and N = 150 for q = 2. The numerical 
results are coincide with that in jUllTj . and the blow-up time is 2.039708648680643 at the first source 
X = 4.079417297361286, corresponding to the maximum value of Umax = 3.16 x 10^. 
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Figure 4: Mesh trajectories and the profiles of u from t = to t = 1.0 for one source case with A^ = 100. 
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Figure 5: Mesh trajectories and the profiles of u for two sources case with iV = 150. 

Example 3 (Sin- type moving sources). We now consider two sources case, in which the sources 
move periodically with the same speed while separated by a constant distance di = 2.5, that is, 

a'o(t) = a'^{t) = Acos{TTt), ao(0) = 0. 

The Blow-up phenomenon is studied in [19j for different amplitudes A. Here we only give the 
numerical results for A = tt (see Figure E]), since all results are similar to those in [TU]. All parameters 
are chosen the same as those in last example. The blow-up occurs at t = 1.689611393639939 on the 
second source with the maximum value of Umax = 3.16 x 10^. 
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Figure 6: Mesh trajectories and the profiles of u for g = 2, A = vr with N — 150. 





Figure 7: Numerical results for symmetric periodic moving sources with N 



Example 4 (Symmetric periodic moving sources). We consider the case for two sources, which 
move periodically and symmetrically. The motion are described by 



ao(t) = j4cos(7rt), ao(0) 



-2.0, 



and ai{t) = —ao{t), with e.g. A 



TT. 



To our best knowledge, there has no theoretical results for multi-sources with different speeds and 
this is the first time numerically investigating the phenomenon for this case. It is shown in Figure 
[7] that blow-up occurs on both sources at t = 2.496881990359248, corresponding to the maximum 
value of ifmax = 3.16 x 10^. 

If local absorbing boundary conditions ()23p are used, the observed domain can be chosen more 
smaller while the results do not be influenced almost. See Figure [5] as an example, where the 
observed domain is set by [ao{t) — 4.0, ai(t) + 4.0], changed as time evolution. Now blow-up occurs 
on both sources at t = 2.496370241342059 with the maximum value of Umax = 3.16 x 10^. 

5 Conclusions 

In this paper, our work focus on the problem of traveling singular sources with different speeds. A 
new moving mesh method in conjunction with a non-overlapping domain decomposition is proposed 
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Figure 8: Numerical results for symmetric periodic moving sources with [xi, x,.] — [cta{t) — 4.0, ai{t) + 4.0]. 

for solving this problems. The whole domain is splitted into q + 1 subdomains by the q sources, 
whose positions are gotten by a predictor-corrector algorithm. Taking the advantages of the domain 
decomposition, the computation of jump [u] is avoided and there are only two different cases dis- 
cussed in the discretization of the physical PDE. Thus, it is easy for the implementation to solve the 
problems with two traveling sources or more. Moreover, the moving mesh method of MMPDEs can 
be applied into each sub-domain respectively. The second-order of the spatial convergence can be 
proved for the new method under a special time marching implementation. The good performance 
of the new method for the blow-up phenomenon is demonstrated through a number of examples 
with two sources. Furthermore, using the new method, we successfully simulate the solutions of 
two sources with different speeds. To our best knowledge, this is the first time investigation for this 
case. The case of three sources or more can be implemented similarly. 
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